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We investigate the effects of anisotropic perturbations in three-dimensional 0(iV)-symmetric 
vector models. In order to assess their relevance for the critical behavior, we determine the 
' renormalization-group dimensions of the anisotropic perturbations associated with the first few 

, spin values of the representations of the 0(N) group, because the lowest spin values give rise to the 

most important effects. In particular, we determine them up to spin 4 for iV = 2, 3, 4, by finite-size 
analyses of Monte Carlo simulations of lattice O(N) models, achieving a significant improvement 
of their accuracy. These results are relevant for several physical systems, such as density-wave sys- 
tems, magnets with cubic symmetry, and multicritical phenomena arising from the competition of 
different order parameters. 
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I. INTRODUCTION AND SUMMARY 



Many continuous phase transitions observed in nature belong to the 0(A) vector universality classes, which are 
characterized by an A-component order parameter with 0(A) symmetry and the symmetry breaking 0(A)— >0(A— 1). 
The superfluid transition in 4 He, the formation of Bose-Einstein condensates, density wave systems, transitions 
in magnets with easy-plane anisotropy, and in superconductors belong to the XY or 0(2) universality class; the 
I ■ Curie transition in isotropic magnets, zero-temperature quantum transitions in two-dimensional antiferromagnets, are 
examples for the Heisenberg or 0(3) universality class; the 0(4) universality class is relevant for the finite-temperature 
transition in two-flavor quantum chromo dynamics, the theory of strong interactions. See, e.g., Refs. 0, 0] for reviews. 

In the absence of external fields, the phase transition of 0(A)-symmetric vector models is driven by only one relevant 
parameter, which is usually associated with the temperature. The corresponding RG dimension is y± = \ jv where 
v is the correlation- length exponent. The leading odd perturbation, which breaks the 0(A) symmetry, is associated 
with the external field h coupled to the order parameter; it has RG dimension y^ — (d + 2 — rj)/2, where rj is the 
exponent controlling the power-law space-dependence of the two-point correlation function of the order parameter at 
criticality. The asymptotic critical power-law behaviors of 0(A)-symmetric vector models have been determined with 
high accuracy. In Table U we report some of the most accurate estimates of the critical exponents v and rj, and of 
the leading and next-to- leading scaling-correction exponents w and W2, which characterize the dominant corrections 
qq , to the universal scaling. 

In this paper we study the effects of anisotropic perturbations breaking the 0(A) symmetry, which cannot be 
t-H \ related to an external vector field coupled to the order parameter, but which are represented by composite operators 
. with more complex transformation properties under the 0(A) group. An interesting question is whether they change 
the critical behavior, or whether they do not affect it so that the symmetry shown by the critical correlations is 
•^h , larger than that of the microscopic model. This issue arises in several physical contexts. Anisotropy in magnetic 
systems may naturally arise due to the cubic structure of the underlying lattice, giving rise to anisotropic interactions 
terms, see, e.g., Ref. [12|]. The relevance of the anisotropic perturbations determines also the nature of the multicritical 
behavior at the meeting point of two transition lines with different 0(ni) and O (n-i ) symmetries , in particular, whether 
the symmetry gets effectively enlarged to 0(ni + 712), see, e.g., Refs. 0, [H, [Tj]. Another interesting issue is the 
critical behavior of secondary order parameters, which are generally represented by powers of the order parameter 
transforming as higher representations of the 0(A) g roup; t heir critical behaviors can be measured in density wave 
systems, such as liquid crystals [H-[I3) see also Refs. [l8l - l2l1 |. 

Let us consider the general problem of the 0(A)-symmetric theory in the presence of an external field h p coupled 
to a perturbation P. Assuming P to be an eigenoperator of the RG transformations, the singular part of the free 
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TABLE I: Some of the most accurate results for the critical exponents of the three-dimensional O(N) vector universality classes 
with ./V = 2, 3, 4, 5. We report estimates of v and r/, and of the leading and next-to-leading scaling correction exponents, 
obtained by lattice techniques (LT) based on Monte Carlo simulations and/or high-temperature expansions, and by quantum 
field theory (FT) techniques such as high-order perturbative expansions. The results without reference have been obtained in 
this paper. A more complete review of results can be found in Ref. @]. 



N 


method 


V 


V 


UJ 


W2 


2 


LT 
FT 


0.6717(1) [3J 
0.6703(15) [4] 


0.0381(2) [3J 
0.0354(25) [4] 


0.785(20) [3J 
0.789(11) [4] 


1.77(7) [5] 


3 


LT 
FT 


0.7112(5) [6] 
0.7117(5) [7] 
0.7116(10) 
0.7073(35) [4] 


0.0375(5) [6] 
0.0378(5) [7] 
0.0378(3) 
0.0355(25) [4] 


0.773 [8] 
0.782(13) [4] 


1.78(11) [5] 


4 


LT 
FT 


0.749(2) [8] 
0.7477(8) [9] 
0.750(2) 
0.741(6) [4] 


0.0365(10) [8] 
0.0360(4) [9] 
0.0360(3) 
0.0350(45) [4] 


0.765 [8] 
0.774(20) [4] 




5 


LT 
FT 


0.779(3) [10] 
0.762(7) [11] 


0.034(1) [10] 
0.034(4) [11] 


0.790(15) [U] 





energy for the reduced temperature t —> and h p — >• can be written as 

Fsi„ g Ht| d V(VI^), 



(i) 



where y p is the RG dimension of h p , and /(x) is a scaling function. Therefore, the RG dimensions of the anisotropic 
external fields quantitatively control their capability to influence or change the asymptotic critical behavior when 
Vp > o. 

In the field-theoretical (FT) framework the 0(A^)-symmetric vector model is represented by the 0(iV)-symmetric 
Landau-Ginzburg- Wilson theory 



H 



d d x 



1 

if' 



+ -u(f 2 ) 2 + h ■ $ 



(2) 



where <I> is an iV-component real field and h an external field. The anisotropic perturbations are conveniently 
classified @, [22[ using irreducible representations of the O(N) internal group, characterized by the spin value I. 
Let us consider the perturbation P m i defined by the power m of the order parameter and the spin representation I of 
the 0(N) group 



where Q dl ' ai is a homogeneous polynomial of degree I that is symmetric and traceless in the I indices: 
Q«($) = $ a 



Qf ($) = $ a $" - jj5 ab $ 2 



<J> 2 



AT- 



Qbfiac _|_ <jjC£a 



Q<f 6cd (<I>) = $ a $ fc $ c $ <i _ — $ 2 (' ( 5 af, $ c $ d + ^ ac $ fe $ d + ,5 a<i (i) b <j) c + (5 bc $ a $ d + <5 M <i> Q <i> c + s cd $ a $ b 

4 v y iV + 4 v 

($ 2 ) 2 (s ab 5 cd + S ac 6 bd + 5 ad 8 bc ) 



(N + 2){N + 4) 



(3) 

(4) 
(5) 

(6) 
(7) 



etc... The classification in terms of spin values is particularly convenient: (i) under the RG flow the operators with 
different spin never mix; (ii) all parameters h^'i' " 1 associated with the components of "° ! have the same RG 
dimension Y m j. On the other hand, operators with different m but with the same I mix under rcnormalization. 
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TABLE II: Estimates of the RG dimensions Yj of the couplings hi associated with the leading anisotropic perturbations Qi 
for the three-dimensional 0(AQ vector universality classes with N = 2,3,4,5. We report results obtained by various methods, 
such as FT perturbative expansions within d — 3 and e-expansion schemes, and lattice techniques, such as high-temperature 
expansions (HT) and finite-size scaling analyses of Monte Carlo simulations (FSS MC). Notice that in the MC estimates of Y4 
reported in Ref. [24| only statistical errors are explicitly given; the authors write that systematic errors are likely of a similar 
size. 



i\ 


method 


i2 (spin Z) 


Y3 (spin 6) 


14 (spin 4) 


2 


Mo -order e expansion 


1.766(6) [11] 


0.90(2) [25] 


-0.114(4) [26] 




FT 6 th -order d = 3 expansion 


1 r 7dd{ 1 fill 
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rbb MC (this paper) 


1.7639(11) 


r\ oai r/nn\ 

0.8915(20) 


—0.108(6) 





FT 5 t]l -order e expansion 


1 ^onVQ 1 ! rill 
1.790(3) [11J 


n n£ f'Q\ ro £1 
U.9o(oJ \Zo\ 


U.UUo(4J |zo| 




FT 6 th -order d = 3 expansion 


1 en/"*^ ri 1 1 
t.sU(oj [1^ 


C\ Cil ( A\ ro Kl 

U.9((4J 


U.Ulo(Dj IzDj 




HT 


1.76(2) [2JJ 








FSS MC 






-0.0007(29) [24] 




FSS MC (this paper) 


1.7906(3) 


0.9616(10) 


0.013(4) 


4 


FT 5 th -order e expansion 


1.813(6) [11] 


1.04(5) [25] 


0.105(6) [26] 




FT 6 th -order d — 3 expansion 


1.82(5) [11] 


1.03(3) [25] 


0.111(4) [26J 




FSS MC 






0.1299(24) [24] 




FSS MC (this paper) 


1.8145(5) 


1.0232(10) 


0.125(5) 


5 


FT 5 th -order e expansion 


1.832(8) [11] 


1.08(4) [25J 


0.198(11) [11] 




FT 6 th -order d — 3 expansion 


1.83(5) [11] 


1.07(2) [25] 


0.189(10) [11] 




FSS MC 






0.23(2) [10] 



The spin-0 operators are already present in the $ 4 Hamiltonian ([2]): the RG dimension of P2 i o is related to the 
correlation length exponent, Y^o = Vt — while the RG dimension of P4.0 (after an appropriate subtraction to 
cancel the mixing with P2,o) gives the leading scaling correction exponent, indeed Y^o = — uj. The spin-1 perturbation 
is related to the external field coupled to the order parameter, thus Yi ; i = yn- 1 Close to four dimensions, thus for 
small e = 4 — d, Y m j < for I > 5, which implies that the only relevant operators have Z < 4. It is reasonable to 
assume that this property holds up to d = 3. Moreover, near four dimensions we can use standard power counting 
to verify that the perturbation with indices m, / mixes with P m > t i, m! < m, but their RG dimensions are significantly 
smaller. In principle, one should also consider terms with derivatives of the field, but again one can show that they 
are all irrelevant or redundant. 

The above arguments show that the most interesting anisotropic perturbations are represented by the spin-2, spin-3 
and spin-4 operators 

Qab r>ab r\abc pa&c r\abcd pabed f a\ 

because they provide the leading effects of anisotropy for each spin sector. As we shall see, the leading RG dimensions 
within each spin sector, 

Yi=Y lyl , (9) 

characterize interesting critical behaviors in various physical contexts. Some Yi have been already estimated by using 
FT approaches based on high-order perturbative calculations, and lattice techniques, such as high-temperature (HT) 
expansions and Monte Carlo (MC) simulations. In Table [TTI we report some results for N = 2,3,4,5. In most cases 
these results provide already a clear indication of the relevance of the perturbation, with the only exception of the 
spin-4 perturbation in the 0(3) universality class, where the value of Y4 is close to zero. While high-order FT results 
indicate the relevance of the spin-4 perturbation, the MC estimate of Y4 appears compatible with zero. Since the 



The perturbation P£ 1 is redundant |2;I . because a Hamiltonian term containing P'i,i can be always eliminated by a redefinition of the 
field 4> a . Anyway, using the equation of motion, one obtains Y3 1 = (d — 2 + rj)/2. 
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issue concerning its relevance is of experimental interest, an accurate determination of I4 is called for to conclusively 
settle it. 

In this paper we present new accurate estimates of the RG dimensions Y; of the anisotropic perturbations for 
N = 2,3,4. For this purpose we perform finite-size scaling (FSS) analyses of Monte Carlo (MC) simulations of 
lattice O(N) spin systems. We achieve a significant improvement of the accuracy of the estimates of Y, essentially by 
combining the FSS method of Rcf . [24] with the use of improved Hamiltonians 28] , which are characterized by the 
fact that the leading correction to scaling is suppressed in the asymptotic expansion of any observable near the critical 
point. Our results are also reported in Table [TT1 As we shall explain later, the errors in the estimates of Yi, and in 
particular of Y4, are quite prudential, they are largely dominated by the systematic error arising from the necessary 
truncation of the Wegner expansions 122] which provide the asymptotic FSS behavior of the quantities considered. The 
results are a good agreement with the estimates obtained by the analyses of high-order FT perturbative expansions, in 
particular with those obtained by resumming 6 th -order d — 3 expansions. Our results show that spin-4 perturbations 
in three-dimensional Heisenberg systems are relevant, with a quite small RG dimension Y4 = 0.013(4), which may give 
rise to very slow crossover effects in systems with small spin-4 anisotropy. The apparent discrepancy with the MC 
result of Ref. [Hj], obtained using the standard nearest-neighbor 0(3) spin model, can be explained by the presence of 
sizable scaling corrections. We overcome this problem by using improved lattice Hamiltonians. The relevance of the 
spin-4 perturbations is important for systems with cubic perturbations |12j |. and also systems whose phase diagram 
presents two transition lines, XY and Ising transition lines, meeting at a multicritical point jl3]. We shall further 
discuss these physical applications later. 

The remainder of the paper is organized as follows. In Sec. [TT| we present the lattice <f) A spin model which we 
consider in our MC simulations, and provide the definitions of the quantities that we consider in our FSS analyses, in 
particular, those related to the spin-Z anisotropics. In Sec. IHII we describe our FSS analyses of MC simulations which 
lead to our final estimates already reported in Table [TTJ Finally, in the conclusive Sec. IIVI we discuss a number of 
physical applications of our results. App. [A] and [B] contain some details of the MC simulations, and further results 
on the critical behavior of O(N) vector models. 

II. THE LATTICE MODEL AND THE ESTIMATORS OF THE ANISOTROPY RG DIMENSIONS 



In this numerical study of O(N) vector models with N — 2,3,4, we consider the </> 4 0(Y)-symmetric lattice 
Haniiltonian 



where <p x is an TV-component real variable, x and y denote sites of the simple-cubic lattice and {xy) is a pair of 
nearest-neighbor sites. In our convention, the Boltzmann factor is given by exp(— Tl^i). For A = we get the 
Gaussian model, while in the limit A — » 00 the 0(Y)-symmetric non-linear a model is recovered. For any < A < 00 
the model undergoes a continuous phase transition in the universality class of the 0(Y)-symmetric vector model. 

In our FSS analyses we consider cubic L 3 lattices with periodic boundary conditions. We consider standard finite- 
volume quantities such as the magnetic susceptibility and second-moment correlation length related to the two-point 
function G{x — y) = (4> x ■ <f> y ), i.e 



A. Improved lattice 0(iV)-symmetric 4> 4 models 




(10) 



(xy) 



x 





(ii) 



X 



and 




(12) 



Another standard quantity for FSS analyses is the quartic Binder cumulant 



Ui 



((M 2 ) 2 ) 
(M 2 ) 2 ' 



(13) 
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The ratio £/L and U4 are RG-invariant phenomenological couplings, thus their large- volume limit at T c is universal. 
We also consider quantities defined keeping one of the phenomenological coupling fixed, in particular keeping the ratio 
£/L fixed, see, e.g., Ref. (2^1 . We define U4 as the Binder cumulant at fixed £/L. 2 

Improved Hamiltonians are characterized by the fact that the leading correction to scaling is eliminated in any 
quantity near the critical point. Therefore in a MC study, the asymptotic behavior at the phase transition can be 
determined more precisely. Improved Hamiltonians were first discussed in Refs. [28[ at the example of the three- 
dimensional Ising universality class using high-temperature series expansions. This idea was first implemented in 
MC simulations of <fi 4 0(7V)-symmetric lattice models for N = 2, 3 and 4 in Refs. @, H3|. In the case of the (f) 4 
lattice model (fTO]) . the improved model is obtained by tuning the parameter A to the particular value A*, where the 
leading 0(L~ U ) scaling corrections vanish in the FSS behavior of any quantity. For this purpose, the RG-invariant 
phenomenological couplings turn out to be particularly useful. Indeed, along the critical line /3 C (A) or keeping another 
phenomenological coupling constant, they behave as 

R{L, A) = R* + c{X)L~ u + ... (14) 

where c(A) is a smooth function of A. Therefore, the equation c(A*) = determines A*. 

The best estimate of A* for N — 2 is A* = 2.15(5) obtained in Ref. 0]. In the case of N = 3, 4, the MC simulations 
performed for this numerical work lead to a revision of the earlier estimates of A*, see App.[B]for details. We obtain 
A* = 5.2(4) for N = 3 and A* = 20±g 5 for N = 4, which update earlier estimates, respectively A* = 4.6(4) of Ref. [|[ 
and A* = 12.5(4.0) of Ref. [8J. 

B. Anisotropy estimators 

In order to compute the spin-Z RG dimensions Y/, we consider appropriate anisotropy correlators. We use the 
magnetization M a — J2 X 4>% ancl the normalized magnetization m a defined as 

M a 

m a = - — r , (15) 
\M\ ' v ' 

to construct objects with given spin properties, such as Q2 b (m), Q3 fcc (m), and Q^^m), obtained by replacing $ a 
with m a in the expressions of Qi, cf. Eqs. ([5]), (O, and ([7]). Then we consider the correlators 



C 2 = ^(^Qf(^)Qf(m)\, (16) 

ab \ x I 

C3 =£(£ Qf c (<t>*)Qf c (m) \ , (17) 

abc \ x I 

C 4 = E ( E Qt cd ^x)Qt cd (m) ) , (18) 



abed \ x 



where Qi(<fr x ) are the operators ((5|), ((6]), and ((7]) constructed using the lattice variable Note that they can be 
rewritten in term of the angle a x defined as <j) x ■ m = \ <j) x cos a x , as 




2 In previous studies, see Refs. [3. [til. [2Sj . another RG-invariant quantity turned out to be very useful, i.e. the ratio Z a /Z p of partition 
functions of a system with anti-periodic boundary conditions in one direction and periodic ones in the other two directions and a system 
with periodic boundary conditions in all directions. Since here we focus on the anisotropy, we have not implemented it to keep the 
project manageable. 
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This expression of C4 shows that it is equal to the improved quantity considered in Ref. [24| to compute the RG 
dimension of the cubic-symmetric perturbation, apart from a constant factor. The asymptotic power-law FSS behavior 
of Ci at T c , i.e. 

Ci ~ L Y > , (19) 

allows us to estimate the RG dimension Yi of the anisotropy associated with Qi . Alternative estimators analogous to 
Ci are also 



. , (V Qf b -(6 x )Qf-(M)) 
D l = Y^ K U2 . Dt~L Y '. (20) 



Note that Q° b -(m) and (Q° b "(iW"))/(M 2 )'/ 2 are by construction RG-invariant quantities (with special symmetry 
properties). Their derivatives with respect to h p , cf. Eq. (JTJ), provide the correlators C; and D\. We also consider the 
corresponding quantities, C/ and D/, at a fixed value of £/L. 



III. FSS ANALYSES OF THE ANISOTROPY CORRELATORS 



In this section we present FSS analyses of high-statistics MC simulations for the 0(2), 0(3) and 0(4) cf> 4 lattice 
models (fT0|) . for values of the parameter A close to A* providing the suppression of the leading scaling correction. 
App. IA II presents some details of the MC algorithm used in the simulations; App. IA 21 reports the values of the 
parameters considered in our MC simulations, the lattice sizes, and the statistics; finally in App. IA3l we discuss the 
behavior of the variance of the observables considered, which influenced the strategy of our FSS analyses of MC 
simulations. 

Most simulations were performed for the 0(3) case, where the spin-4 RG dimension I4 is close to zero, and therefore 
high accuracy is needed to determine its sign. This task is made particularly hard by the rapid increase of the cost to 
get accurate data for C4 and D4 with increasing the lattice size, essentially due to a significant increase of their variance, 
see the discussion in App. IA 31 As a consequence, our FSS analyses to determine I4 are limited to relatively small 
lattice sizes. On the other hand, the systematic error due to the necessary truncation of the Wegner expansion [22| . 
see Eq. (|2"Tj) below, of the quantities considered turns out to be significant, and its reduction requires accurate results 
for large lattice sizes. This represents the major limitation for the accuracy of our numerical determination of Y4. 

App. [B]reports further FSS analyses of the MC simulations which allow us to update some of the results concerning 
the 0(N) vector models, such as the estimates of A*, of the critical exponents and other universal quantities. 



A. General strategy of the FSS analysis 

In order to obtain accurate estimates of the universal quantities, such as the critical exponents and RG dimensions 
Yi, it is important to have a robust control of the corrections to the asymptotic behaviors, which are suppressed by 
powers of the lattice size L. The behavior of general quantities introduced to estimate critical exp onents, such as C; 
and Di defined in the previous section, can be expressed by an asymptotic Wegner expansion [22] as 

A(X; L) = c(X)L y [l + a(A)L~ w + ^ a^L'^] (21) 

i=2 

where y is the leading universal exponent that one wants to accurately estimate. In the case of 0(2), 0(3) and 0(4) 
vector models the leading scaling correction exponent is given by w k 0.8, see Table [U Numerical approaches based 
on improved Hamiltonians allow us to suppress these leading scaling corrections, and also those related to nu, where 
n = 2, 3, 4, whose coefficients behave as (A — A*)". The next-to-leading correction is controlled by the exponent cj 2 , 
estimated in Ref. [5( by L02 ~ 1.8, see Table HI Then there are well established corrections with uji ~ 2, for example 
related to the breaking of spatial rotational invariance in cubic lattice systems [3l| . but also to analytic backgrounds, 
etc... Moreover, in the case of the spin-Z anisotropy correlators, we may also have scaling corrections induced by 
higher-dimensional spin-Z operators, such as P; + 2./, cf. Eq. ([3]). On the basis of a dimensional analysis around four 
dimensions, they are expected to give rise to scaling corrections suppressed by powers k; = 2 + 0(e), as also shown 
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by the 0(e) calculation of the difference of the RG dimensions of the anisotropy operators Pi+2.1 and Pij, which is 3 

Yi + 2,i - Y u = -2- e6(l - 1)/{N + 8) + 0(e 2 ). (22) 

In known cases for the spin-0,1,2 sectors, the difference between RG dimensions of the same sector remains close 
to their four dimensional values. Therefore, as a prudential procedure, after curing the residual 0(L _CJ ) scaling 
corrections, see also below, we must consider possible 0(L~ K ) scaling corrections with k > 1.6. 



1. Residual leading scaling corrections in approximately improved Hamiltonians 



Residual leading scaling corrections are generally present due to the fact that A* is only known approximately, and 
also because the MC simulations are usually performed close but not exactly at the best estimate of A*, which is 
usually determined at the end of the MC simulations. For example, in the case N = 3 our best estimate is A* = 5.2(4), 
while most MC simulations were performed at A = 4.5, and others at A = 4 and A = 5 for smaller lattices to determine 
A*. 

The residual 0(L _W ) corrections, due to the fact that A is close but does not coincide with its optimal value A*, 
can be further suppressed as follows. The basic idea is that leading corrections to scaling can be best detected by 
analyzing the Binder cumulant U4 at a fixed value of £/L. At a generic A = Ao we have 

U4X ;L) = U:+au(X )L- u + ..., (23) 

where U\ is the universal large- volume limit on a periodic L 3 box at fixed £/L, which of course depends on which 
value of £/L is chosen. Then, we consider a pair Ai, A2, where one of the two values may be equal to Ao, and the 
differences 



A l/ (A 1 ,A 2 ;L) = C7 4 (A 2 ;L) - C7 4 (A i; L) 
where the leading large-volume contributions cancel, thus they behave as 

A c/ (A 1 ,A 2 ;L) = 6 c/ (Ai,A 2 )L- w + ... . 
The amplitude 6[/(Ai, A 2 ) = aj/(A 2 ) — ajj(Xi) can be estimated by fitting the data to 

A(X 2 ;L) 



(24) 



(25) 



r A {Xi,X 2 ;L) 



A(X i; L) 



Finally, we take ratios 

(26) 



of the quantity A that we intend to correct to eliminate the residual 0(L u ) corrections. Their data can be fitted to 
its large- L behavior 



r A (Xi, A 2 ; L) = ^\ [l + 6(Ai, X 2 )L~ 
c(Ai) 



(27) 



where 6(Ai,A 2 ) = a(A 2 ) — a(Ai) and a(X) is the amplitude of the 0(L _ ") corrections, cf. Eq. (f2"Tj). Notice that 
it is simpler to extract 6(Ai,A 2 ) than a(A) from the numerical data, because, beside the cancellation of the power 
divergence L v , also subleading corrections cancel to a large extent. Now we use the universality of ratios of correction 
amplitudes, which implies 



a(A ) _ b(X 1 ,X 2 ) 



au{Xo) bu(Xx,X 2 ) 
In order to eliminate the leading 0(i _w ) corrections from A, we construct 

6(Ai,A a ) 



1 A (X 0] L) = A(X ;L) 



1 - 



bu(Xi,X 2 ) 



au(X )L~ 



(28) 



(29) 



We note that within e expansion the operator P;+2 ; mixes with other spin-i operators containing derivatives (two derivatives instead 
of $ 2 ), but this mixing contributes to 0(e 2 ). 
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A lnC 2 
□ lnC 5 




1 2 3 4 5 6 

InL 

FIG. 1: (Color online) Log-Log plots of C*2 and C3 versus L at /3 C and A = 4.5. The errors of the data are hardly visible. 

This procedure eliminates the leading 0{L~ U ) scaling corrections, allowing us to neglect them in the fits of the data 
of Xa(Ao; L) to estimate the leading exponent y. 4 

We also mention that alternative procedures, based on the idea of defining improved observables with suppressed 
leading scaling corrections, are outlined in Refs. pit l32|. 

2. Next-to -leading corrections 

Next-to- leading corrections arise from the term associated with UJ2 ~ 1.8, and the others with exponents close to 
two. In the fits of the data, even with high statistics data as we have here, only a very limited number of correction 
terms can be taken into account. The truncation of Eq. (|21[) leads to systematic errors in the results for the exponent 

V- 

One way to control these systematic errors is to study several quantities A^ that have the same critical behavior: 

A( n \L) = c n Ly(l + J2amL- Ui ) (30) 

i 

In general one might expect that for different A' 71 ' the coefficients a n i are different. Therefore the variation of the 
estimate for y obtained by fitting several provides an estimate of the systematic error. However, in our case we 
have only the two quantities Ci and Di, which are closely related. Therefore we would like to estimate the systematic 
error by fitting a single quantity. To this end we consider the Ansatz 

A(L) = cL y (l + aL' u + a 2 , c ffi^ 2 ' otf ) (31) 

(for improved models a = 0), with 

w 2 ,cff > 1.6 (32) 

Barring an unlike significant cancellation between different correction terms, there must be a value of W2,cff > 1.6 such 
that y takes its correct value. Since we expect that, as long as correction are small, the resulting y is a monotonic 
function of k>2,eff) we use the results obtained for u>2,eS = 1-6 and W2, e ff = 00 (i.e. without the term C2 je ff£ - " 2,otf ) as 
bounds for the correct result for y. 



4 The coefficient c = ajy (Ao)fc(Ai , A2 )/!>[/ (Ai , A2) is numerically determined with an error Ac, which is usually dominated by the uncertainty 
on d[/(Ao). This error can be taken into account by computing Xa(Ao; L) using c and c ± Ac. The difference between the results of 
their fits is essentially related to the error due to the uncertainty of our estimate for A* , since also the uncertainty of the estimate of A* 
is mainly caused by the error of aj/(Ao), see also App. IB 21 
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FIG. 2: (Color online) Log-Log plots of C4 versus L at /3 C and A = 4.5 



B. Results for the spin-i RG dimensions 

1. The 0(3) model 

To begin with we present the FSS analysis of the data for the 0(3) model. In order to give an idea of the quality 
of our data, we show the data of C; at /3 C and A = 4.5 in Figs. [T]and[2J In all cases, including C4, the data clearly 
increase with increasing L, indicating the relevance of the perturbation. Note that the error of C4 is rapidly increasing 
with increasing L, see App. IA3l for details. 

We analyze various quantities to estimate the RG dimensions Yf. the original quantities C; and Di introduced 
in Sec. Ill Bt their counterpart Ci and Di computed at the fixed value £/L = 0.5644 (which is a good estimate of 
the large-volume limit of £/L at j3 c , see App. |B|, and also the quantities Ci^ mp — UfCi and A, imp = UfDi a gain 
taken at £/£ = 0.5644 where the exponent x is chosen to further suppress the leading corrections (see Refs. [1,|32[ 
for details). In principle, the latter quantities should be more suitable for the numerical analysis. Indeed, by fixing 
£/L = 0.5644 we avoid the error due to the uncertainty of j3 c , and by the construction of the improved observables the 
effect of the uncertainty of A* is strongly reduced. However also subleading corrections vary, and, unfortunately, they 
become numerically larger in these cases. Nethertheless it is useful to study these quantities. Since the amplitudes 
of corrections change, these modified quantities give us additional control over the systematic error that is caused by 
truncated ansaetze. 

Let us now discuss the analysis of the quantities Di in some detail. In the case of the quantities C; we proceed in 
a similar way. Following the discussion of section IIII A II we first analyze the ratios 

Di{\ = 5, fi = 0.687564) 
r ° l A(A = 4,0 = 0.68439) 1 ' 

where j3 — 0.687564 and (3 — 0.68439 are the estimates for (3 C given in Table V of Ref. 6] and 

= A(A = 5) 
r °< A(A = 4) ■ 

We fit these ratios to the ansatz 

r = c(l + Wr w ) (35) 

where we set uj = 0.79. To give an idea how accurately the coefficient b can be determined let us discuss a few 
examples. In the case of D2 a fit of all data with L > L m i n with L m i n = 6 gives the result b — —0.00841(38) 
and x 2 /DOF= 3.80/9. Increasing L m i n the estimate of b changes very little, for example, for L m i n = 8 we obtain 
b = -0.00830(64) and x 2 / D F = 3.58/7. In the following analysis we shall assume b = -0.0083(7). In the case of D 2 
we obtain for L mm = 7 the result b = 0.00348(24) and x 2 /DOF= 8.48/8 and for L mm = 9 the results b = 0.00407(41) 



(34) 
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TABLE III: Fits of Xn 2 (column 2 and 3), Iq 2 (column 4 and 5) and U4D2 (column 6 and 7) with the ansaetze (J36J) and (137) . 
We give the L m in of the fit, which is typically the smallest L m i n that produces an acceptable fit and the result for Y<z. 



ansatz 


Lmin 


Y 2 


Lmin 


Y 2 


Lmin 


Y 2 


p6|) 


24 


1.79067(5) 


28 


1.79078(3) 


32 


1.79080(5) 


(57) 


12 


1.79019(7) 


8 


1.79053(2) 


8 


1.79049(2) 



TABLE IV: Fits of Td 4 (column 2 and 3), If, 4 (column 4 and 5) and £7f D4 (column 6 and 7) with the ansaetze (J36J) and l[3"?]l. 
We give the L m i n of the fit, which is typically the smallest L m i n that produces an acceptable fit and the result for Y4. 



ansatz 


Lmin 


Y 4 


Lmin 


Y A 


Lmin 


Y 4 


([36]) 


14 


0.0143(8) 


14 


0.0142(8) 


16 


0.0160(10) 


m 


12 


0.0122(26) 


12 


0.0127(25) 


12 


0.0122(26) 



and x 2 /DOF= 2.34/6. In the following analysis we shall assume b — 0.004(1). It is interesting to observe that, 
by taking D2 at £/L = 0.5644 instead of /3 C , even the sign of the correction amplitude changes. For D4 we obtain 
b = 0.0313(45) and x 2 /DOF= 7.49/8 using L m i n = 7. The result changes little when we increase L m i n . For example, 
we get b = 0.0303(93) and % 2 /DOF= 7.42/6 for L min = 9. In the following we shall assume b = 0.03(1). For D 4 we 
get instead b = 0.05(1). Note that also the correction amplitudes of D4 and D4 are different. 

In order to compute the quantities X]j l and Xq , defined as in Eq. (|29[) to suppress the residual leading scaling 
corrections, we use 6^/(5,4) = —0.01126(4) and ay (4.5) = 0.007(4) as obtained in appendix IB 21 In the product UfDi 
the choice x = —U^b/bu eliminates leading corrections to scaling. The advantage of this quantity is that it does not 
require au, which is affected by a relatively large error. 

Next we have fitted the resulting quantities with the ansaetze 



1 MAi,A 2 ) 



= aL Yl (36) 



l Dl (\o;L) = A(A ;i) 

and 

X Dl =aL Y '(l + dL- 1 - 6 ) (37) 

and correspondingly for the quantities and t/f A- The effect of the uncertainties of /3 C , and the quantities au, bu, 
b\ need to construct Xp^ X^ i and U4D1, are estimated by varying these input parameters. E.g. in order to estimate 
the uncertainty of Xq induced by the uncertainty of au, we have repeated the fits using data where we have used in 
eq. (|29p the central value of au plus its error instead of the central value. 

In Table Hill we report results of fits for Xd 2 , Xq 2 , and t/J-C^. We note that the estimates of Y2 obtained by the 
two fits and the three quantities differ by larger amounts than their statistical errors. Hence systematic errors are 
more important than the statistical one. Taking into account also the results obtained for C4 and the quantities 
derived from it we arrive at our final estimate Y2 — 1.7906(3) which covers most of the acceptable fits and also takes 
into account the uncertainties in the construction of X]j 2 , Xq and Tj%Di. In a similar way we arrive at the estimate 
Y 3 = 0.9616(10) of the spin-3 RG dimension. 

Finally, let us discuss the analysis leading to our estimate of Y 4 . In Table HVl we give some results of the fits with 
the ansaetze (136)) and (l37l) . As our final result we quote Y 4 = 0.013(4) which covers all estimates given in Table HVl 
The uncertainties in the construction of Xd 4 , Xf )i and U4D4 are taken into account. Furthermore, this estimate is 
fully consistent with the results obtained from the analysis of Xc 4 , Xq a and U4C4. 

We conclude with a few remarks on the possibility of further improving the estimate of Y4. Its accuracy is essentially 
limited by the fact that the variances of the correlators C4 and D4 rapidly increase with increasing lattice size, not 
allowing us to get accurate results for large lattices, indeed extremely high statistics are necessary for L > 32 already. 
Thus, the reduction of the systematic error due to the truncation of the Wegner expansion appears quite problematic, 
because it can only get reduced by accurate results for larger lattice sizes. One purely technical way in this direction 
could be the simulation with local algorithms (Metropolis + many overrelaxation sweeps) on GPUs (Graphics cards). 
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2. The 0(2) and 0(4) models 

In the cases of the XY and 0(4) universality classes we have determined the exponents along similar lines, obtaining 
the results reported in Table [TT We only mention that, since in the case of the XY universality class, A* and fa at 
A = 2.1 are accurately known 3|, we abstained from analyzing the quantities U%C\ and U%Di. In the case of the 
0(4) universality class the situation is different; here we do not have a very precise estimate of A* and also fa is 
only moderately well known at A = 12.5, where most of our simulations are performed. Therefore we have based our 
analysis on Ci and Di and the improved quantities UfCi and UfDi, where the quantities are taken at £/L = 0.547. 



IV. CONCLUSIONS AND DISCUSSION OF SOME APPLICATIONS 



In this paper we study the effects of anisotropic perturbations in three-dimensional 0(iV)-symmetric vector models, 
which cannot be related to an external vector field coupled to the order parameter, but are represented by composite 
operators with more complex transformation properties under the O(iV) group. For the models with N = 2,3,4, 
we determine the RG dimensions Yj of the anisotropic perturbations associated with the first few spin values of the 
representations of the O(N) group, because the lowest spin values give rise to the most important effects. This is the 
first numerical study based on MC simulations for the spin-2 and spin-3 perturbations, while MC results for spin-4 
operators were already reported in Ref. [24| . 

We present FSS analyses of MC simulations of improved Hamiltonians with suppressed leading corrections to 
scaling, which allows us to achieve a robust control of the systematic errors arising from scaling corrections. Our 
results are reported in Table HH together with earlier results by various approaches. They are in good agreement with 
the estimates obtained by field-theoretical methods, by resumming high-order perturbative series. Our results show 
that spin-4 perturbations in three-dimensional Heisenberg systems are relevant, with a quite small RG dimension 
I4 = 0.013(4), which may give rise to very slow crossover effects in systems with small spin-4 anisotropy. 

In the following we discuss a number of physical systems where the results of this paper for the anisotropic pertur- 
bations can be used to infer the critical behavior of some physically interesting quantities. 



A. Critical exponents of secondary order parameters 

Beside the standard critical exponents associated with the order parameter, density wave XY systems allow to 
measure the higher- harmonic critical exponents related to secondary order parameters, which can be theoretically 
represented by polynomials of the order parameter with spin representation higher than one, such as the spin-Z 
operators Qi{4> x ), cf. Eqs. ([S][7]). 

The behavior at zero- momentum of the correlation functions involving the operators Qi{4> x ) can be described by 
introducing an appropriate external field hi coupled with Qi(<f> x ), and writing the singular part of the free energy as 
in Eq. (TTJ). Then, differentiating with respect to hi, we obtain the behavior of the secondary magnetizations in the 
broken phase, 

(Ql(<t>x)) ~ \tf l , Pi=v{d-Yi). (38) 

Our estimates of the RG dimensions Yi for the XY universality class, Y2 = 1.7639(11), Y3 = 0.8915(20) and Y4 = 
-0.108(6), give 

fa = 0.8303(8), fa = 1.4163(13), fa = 2.09(4). (39) 
Moreover, the nonanalytic scaling behaviors of spin-Z susceptibilities are 

Xi =^(Q i (0o)gK0x))^|ir 7! , 71 = "(2*1 - d), (40) 

with 

72 = 0.3545(15), 73 = -0.817(3), 74 = -2.160(8). (41) 

Note that the power law \t\~ 11 in the susceptibility \i represents the leading term only if 7; > 0, otherwise the 
nonuniversal analytic contributions provide the dominant behavior, see, e.g., Ref. [HJ]. We also mention that the 
structure factor, obtained by Fourier transforming the correlation function Gi(x — y) = (Qi{4> x )Qi(4> y )} , is expected 
to behave as Gi(q) ~ |i|~ 7i /;(?£)> where f\ is a universal function, see Ref. [33[ and references therein. 
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Discussions of the experimental systems and results for the higher-harmonic exponents can be found in Refs. [2], 
[25|, The experimental estimates are in substantial agreement with the theoretical results. Here we only mention 
a few of them. Analyses [HI, [l(| S] of the experimental data near the smectic-C-tilted-hexatic-I transition provided 
estimates of the crossover exponent <fri — Yiv. By replacing v = 0.6717, they give I2 = 1-7(1) and Y3 = 0.6(3). In 
Rcf. 20] the estimates (3 2 = 0.87(1) and (3 3 = 1.50(4) were obtained for Rb 2 ZnCl 4 . 



B. Magnets with cubic symmetry 



The magnetic interactions in crystalline solids with cubic symmetry, like iron or nickel, are usually modeled by 
using the 0(3)-symmetric Heisenberg Hamiltonian with short-range spin interactions, such as 

•ffspin = — J E Si ■ Sj (42) 

where S 2 — 1 and the sum is over nearest neighbors. However, this is a simplified model, since other interactions 
are present. Among them, the magnetic anisotropy that is induced by the lattice structure (the so-called crystal 
field) is particularly relevant experimentally, see, e.g., Ref. [351 ] . In cubic-symmetric lattices it gives rise to additional 
single-ion contributions, the simplest one being 

EE 5 ? 4 - ( 43 ) 

i a 

These terms are usually not considered when the critical behavior of cubic magnets is discussed. However, this is 
strictly justified only if these nonrotationally invariant interactions, that have the reduced symmetry of the lattice, 
are irrelevant in the RG sense. The corresponding cubic-symmetric perturbation X^a^" 4 *° the O(iV) theory is a 
particular combination of spin-4 operators P% b £ d and of the spin-0 term P-^o, 

N 

E $a4 = E pvrm + ^^ot*) (44) 

a a— 1 

Since i-4,0 is always irrelevant, the relevance of the cubic-symmetric anisotropy is related to the value of the spin-4 
RG dimension Y4, and in particular to its sign. Our results, and in particular Y4 = 0.013(4) for the 0(3) universality 
class, show that the cubic perturbation is relevant at the three-dimensional O(N) fixed point when N > 3, confirming 
earlier FT results [26|, I36l438j . This implies that for N > 3 the asymptotic critical behavior is described by another 
cubic-symmetric fixed point, see, e.g., Refs. @ for a general discussion of the RG flow in the $ 4 theories with cubic- 
symmetric anisotropy. However, differences between the Heisenberg and cubic critical exponents are very small [llj ]. 
for example v differs by less than 0.1%, which is much smaller than the typical experimental error for Heisenberg 
systems !2j. Therefore, distinguishing the cubic and the Heisenberg universality class is very hard in experiments. 



C. Multicritical phenomena in 0(ni)®0(ri2)-symmetric systems 



The competition of distinct types of ordering gives rise to multicritical behaviors. The multicritical behavior arising 
from the competition of two types of ordering characterized by O(n) symmetries is determined by the RG flow of the 
most general 0(ni)©0(n 2 )-symmetric LGW Hamiltonian involving two fields (f>\ and </> 2 with n\ and n 2 components 
respectively, i.e. [13| 



rl 



1 



1 



1 



(45) 



A multicritical point (MCP) is achieved when r\ and r 2 are tuned to their critical value, and the corresponding 
multicritical behavior is determined by the stable FP of the RG flow of the quartic parameters. It may occur at the 
intersection of two critical lines characterized by different O(ni) and 0(n 2 ) order parameters. 

An interesting possibility is that the stable FP has 0(ni + n 2 ) symmetry, so that the symmetry gets effectively 
enlarged approaching the MCP. The stability properties of the 0(ni +n 2 ) symmetric FP can be inferred by noting [ll| 
that the Hamiltonian (|45p contains combinations of spin-2 and spin-4 polynomial operators with respect to the 
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0(ni + 712) group, which are invariant under the symmetry 0(ni)©0(ri2). Defining $ as the [nx + n2)-component 
field </> 2 ), they are given by the spin-0 operators "J 2 and ( ( & 2 ) 2 , by the spin-2 operators 



ni 



O2.2 = V P$% =<t>t ^* 2 , 4 ,2 = $ 2 2 , 2 , (46) 



0=1 



and by the spin-4 operator 



04,4-^ ^ p m -M,- ni + n2 + 4 + (ni + n2 + 2)(ni + n2 + 4) - ( 4? ) 



a— 1 + 1 



The 0(ni + 712) FP controls the multicritical behavior if it is stable against the fourth-order perturbations, and, 
in particular, the dominating spin-4 perturbation 04,4, (the perturbation O42 is expected to be irrelevant after the 
subtraction of its lower-dimension spin-2 content [ll(). 

Our FSS MC results for the spin-4 RG dimensions Y4 (see Table HI]) . and, in particular, that for the 0(3) uni- 
versality class, provide a conclusive evidence that Y4 > for n\ + 112 > 3, confirming earlier indications from FT 
computations [ll|. Therefore the enlargement of the symmetry 0(m)©0(n2) to 0(ni + 0,2) does not occur, unless 
an additional parameter is tuned beside those associated with the quadratic perturbations. We may observe an en- 
largement of the symmetry to 0(2) only when two Ising lines meet. In this case the RG dimension Y2 of the spin-2 
operator 02,2 provides the crossover exponent cf) = 1/Y2 = 1.1848(8) at the MCP. 

These results can be applied to the study of the phase diagram of anisotropic antifcrromagnets in a uniform magnetic 
field H\\ parallel to the anisotropy axis, which present a MCP in the T — H\\ phase diagram, where two critical lines 
belonging to the XY and Ising universality classes meet [HI, Il4j |. Experimental realizations of these systems are 
reported in Refs. [39l - l4lT |. which typically show phase diagrams with a bicritical MCP. The initial hypothesis of an 
enlarged 0(3) symmetry at the MCP, on the basis of low-order FT calculations [l4|], was then questioned by high- 
order FT computations [ll| (see also Ref. [42j]), indicating a very weak instability of the 0(3) FP. This instability 
was then questioned by the numerical MC study of Ref. [43j |. where evidence of a 0(3)-symmetric bicritical point is 
claimed in the phase diagram of the so-called XXZ model, which models anisotropic antiferromagnets in an external 
field, showing a MCP where an XY and an Ising transition line meet. Actually, this result was one of the major 
motivation of this numerical work to further check the relevance of the spin-4 perturbation at the 0(3) FP, because 
an asymptotic 0(3) multicritical behavior requires I4 < 0. Our MC results fully confirm earlier high-order FT results, 
i.e. the relevance of the spin-4 0(3)-breaking term which are generally present in these models. This implies that 
a bicritical point in the Heisenberg universality class is excluded, unless one achieves a complete cancellation of the 
spin-4 term by an appropriate fine tuning. 

As inferred by FT calculations, the actual stable FP has a biconical structure [ll|. A quantitative analysis of the 
biconical FP shows that its critical exponents are very close to the Heisenberg ones. For instance, the correlation- 
length exponent v differs by less than 0.001 in the two cases. Thus, it should be very hard to distinguish the biconical 
from the 0(3) critical behavior in experiments or numerical works based on Monte Carlo simulations. 

The crossover exponent describing the crossover from the unstable 0(3) critical behavior is very small, i.e. 04 = 
vY± = 0.009(3), so that systems with a small effective breaking of the 0(3) symmetry show a very slow crossover 
towards the biconical critical behavior or, if the system is outside the attraction domain of the biconical FP, towards 
a first-order transition. Thus, they may show the eventual asymptotic behavior only for very small values of the 
reduced temperature. Likely, the numerical analysis of Ref. [43| was just observing crossover effects. 
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Appendix A: Monte Carlo simulations 

1. Monte Carlo algorithm 

As Monte Carlo algorithm we use a hybrid of the local Metropolis, the local overrelaxation and the single cluster 
[44| algorithm. The proposals for the local Metropolis update are given by 

(f>' x = <f> x + sr x (Al) 
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where s controls the step size and the components of the random vector r x are uniformly distributed in the interval 
[—0.5,0.5]. This proposal is accepted with the standard acceptance probability 

P acc = min[l,exp(-AH)] . (A2) 

The step size s is chosen such that the acceptance rate is roughly 50%. In the case of the local overrelaxation update, 
the new value of the field is given by 

where $ x — ^2 y nn x <fi y is the sum over all fields that live on sites y that are nearest neighbors of x. In the case of the 
local updates we run through the lattice in typewriter fashion. Going through the lattice once is called one sweep. 
We use the following cycle of updates: One Metropolis sweep, one overrelaxation sweep, L/2 single cluster updates, 
two overrelaxation sweeps and finally L/2 single cluster updates. In this cycle, we compute the observables after L/2 
single cluster updates, i.e. twice. 

The average size of a cluster is proportional to the magnetic susceptibility that grows like L 2 ~ v . Therefore, with 
our choice of L/2 single cluster updates per cycle, the fraction of sites that is updated by the cluster algorithm in one 
cycle of the algorithm stays roughly constant. We also note that the overrelaxation update takes very little CPU time 
compared with the Metropolis update. For L = 32 and N — 3 the CPU time needed for one overrelaxation sweep, 
one Metropolis sweep, and L/2 single cluster updates roughly behave as 1 : 4 : 3. 

In all our simulations we have used the SIMD-oriented Fast Mersenne Twister algorithm [45| as pseudo-random 
number generator. 

2. Statistics of the simulation 

In the case of the XY universality class, we performed most of our simulations at A = 2.1 and j3 — 0.5091503. We 
simulated the lattice sizes L = 6,7, 8, ...,18 and 20,22,24,26,28. Throughout we performed 10 9 measurements. In 
total these simulations took about 7 month of CPU time on a single core of a Quad-Core AMD Opteron(tm) Processor 
2378 running at 2.4 GHz. In addition we performed simulations at A = 2.2 and j3 = 0.508336 where we simulated the 
lattice sizes L — 6, 7, 8, 12. The results for A = 2.2 are used to estimate the effect of the uncertainty of A*. Note 
that A* = 2.15(5) [J]. The values of /3 chosen for the simulations at A = 2.1 and 2.2 are the estimates of /3 C given in 
Table II of Ref. [|. 

In the 0(3) case we performed most simulations for A = 4.5 which is close to our old estimate A* = 4.6(4) Q. We 
simulated at j3 = 0.686238 which is close to the estimate f3 c = 0.6862385(20) @. For the lattice sizes L = 6,7, 8, 9, 16 
we performed 10 9 measurements, for L = 17, 18, ...,32 between 1.1 x 10 9 and 1.2 x 10 9 measurements and 5 x 10 8 , 
2.5 x 10 s , and 10 6 measurements for L = 48, 64 and 256, respectively. In total these simulations took about 4 years 
of CPU time on a single core of a Quad-Core AMD Opteron(tm) Processor 2378 running at 2.4 GHz. In addition, we 
performed MC simulations at A = 4.0, /3 = 0.68439 and A = 5.0, /3 = 0.687564 on lattices of the size L — 6,7, 8, 16. 
Throughout we performed 10 9 measurements. These results are used to determine our new estimate of A* and the 
effect of the uncertainty of A* on our estimates of the RG exponents. 

In the 0(4) case most of our simulations were done for A = 12.5 and (3 = 0.9095167. For L = 6,7, 8,. ..,18 
and 20,22,24,26,28 we performed 10 9 measurements and for L — 40 we performed 6.5 x 10 s measurements. For 
L = 256 we performed 10 6 measurements and simulated at (3 = 0.909513, which was our preliminary value of j3 c . 
This simulation was done to get a better estimate of (3 C . In this simulation we did not measure the quantities C\ and 
Di. The estimate (3 = 0.9095167 used above was obtained by requiring that £/L = 0.547 which is the result for the 
large volume limit (£/£)* of [8|. In addition, in order to determine A* and the effect of the uncertainty of A* on the 
accuracy of our estimates of the RG-exponents, we have simulated at A = 14 the lattice sizes L = 6, 7, 8, ...12; A = 18 
the lattice sizes L = 6, 7, 8, ...12; A = 22 the lattice sizes L = 6, 7, 8, ...16, 18, 20; A = 30 and 32 the lattice size L = 6; 
and for A = oo the lattice sizes L = 6, 7, 8, ...12, 16, 24, 32. Throughout the statistics is 10 9 measurements. 

The CPU time used for the whole study amounts to roughly 7 years on a single core of a Quad-Core AMD 
Opteron(tm) Processor 2378 running at 2.4 GHz. 

3. Variance of the observables 



The behavior of the variance of the quantities considered in our MC simulations strongly affects the design of our 
study. The main problem, as already observed in ref. 24] is that the relative statistical error, at a fixed number of 
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updates, of C4 and D4 rapidly increases with the lattice size. Therefore we have to focus on smaller lattice sizes than 
one would do in a study mainly aiming at the exponents v and 77. 

Let us discuss this problem in a bit more detail at the example of the simulations for N = 3, A — 4.5 and the 
quantities Di. Since we average over 10000 measurements at simulation time, we can not disentangle integrated 
autocorrelation time and variance of the quantities. Therefore in the following we discuss the relative statistical error, 
normalized to 10 9 measurements. In the case of D4 this relative statistical error is increasing from 0.000175 for L = 6 
up to 0.051 for L = 256. This increase is well described by a power law e oc L x , with x » 1.45. Also in the case 
of Z?3 the relative error is increasing; 0.000064 for L — 6 up to 0.00022 for L — 256. However here the increase is 
smaller; it is characterized by the exponent x » 0.3. Interestingly, for D2 we find that the relative statistical error is 
even decreasing a bit; 0.000037 for L = 6 down to 0.00003 for L — 256. The corresponding exponent is x « —0.05. 
This behavior can be compared with that of the relative error of the slope of the Binder cumulant or the second 
moment correlation length. These quantities are used to determine the critical exponent v. In both cases we find a 
mild increase of the relative error, which is characterized by the exponents x w 0.06 and x ks 0.14, respectively. 

As shown in Ref. [HI, the problem of the large variance of C4 can be reduced by performing a larger number of 
overrelaxation updates which are relatively cheap in terms of CPU time and measure C4 after each such update. This 
way one could improve the efficiency in terms of 1/ [(CPU-time) x error 2 ] of C4 or D4 by about a factor of 2 compared 
with the update cycle used in our simulations. However, since this would have an adverse effect with respect to all 
other quantities that we have measured we abstained from this. 

For several observables, such as the susceptibility and the quartic Binder cumulant, the statistical errors at fixed 
£/L arc smaller than those at fixed j3 close to (3 C . Some comparisons are reported in Refs. [!, l4q |. This is due to 
cross correlations and to a reduction of the effective autocorrelation times. Taking C; or _D; at £/L fixed reduces the 
variance in a /-dependent way. For the C4 and D4 cases there is virtually no reduction of the error. For L = 6 there 
is still an improvement by a few percent, however with increasing L, the ratio of errors goes rapidly to 1. In the I = 3 
case we observe a mild improvement by fixing £/£. For C3 the ratio of statistical errors is 1.9 for L = 6, 1.10 for 
L = 64 and 1.017 for L = 256. In the case of D$, the ratio of statistical errors is 1.33 for L = 6, 1.06 for L = 64 and 
1.014 for L = 256. The reduction of the statistical error is most significant in the I — 2 case. For C2 the ratio of the 
statistical errors is 3.49 for L = 6, it decreases to 2.65 at L = 27 and then increases again; 2.69 at L — 64 and 2.88 
for L = 256. For D2 the ratio of the statistical errors is 2.21 for L = 6, has its minimum 1.91 at L = 23, takes 2.02 
for L = 64 and 2.20 for L = 256. 



Appendix B: Some further results for the O(N) vector models, N = 3 and 4 

1. New estimate for /3 C 

In order to determine /3 C , we fit the data for £/L and U4 at A = 4.5 to the ansaetze 

R(L,f3 c )=R* (Bl) 



R{L, f3 c ) = R* + aL~° (B2) 

and 

R(L, p c ) = R* + aL- 79 + bL~ e (B3) 

where either e = 1.6 or e = 2. Here we take 0.79 as value of the correction exponent u. By replacing it with 0.77 say, 
our results for /3 C and R* change only very little. In this study, we only calculate first derivatives of the quantities; 
therefore in the fits we use the approximation 

R(L,/3)ttR(L,/3 s )+a(p-l3 s ) (B4) 

where f3 s is the value of the inverse temperature used for the simulation. Since j3 s is very close to our final result for 
/3 C , the error due to the truncation of the Taylor-series can be ignored. 

Let us first discuss the analysis of £/L. Taking no corrections into account, i.e. fitting with the ansatz (|B1[) . x 2 /DOF 
remains unacceptably large until most of our lattice sizes are discarded. Including L = 48, 64 and 256, we obtain 
(£/L)* = 0.56421(5), (3 c -f3 s = -0.0000006(5) and x 2 / D OF= 1.72/1. Using the ansatz (|B2j). i.e. adding a correction 
term aL~ - 79 we get a x 2 /DOF smaller than 1 starting from L m i n — 12, where all lattice sizes L > L m i n are taken into 
account. Discarding further data points x 2 /DOF is further decreasing and (£/£)* and /3 C — f) a move monotonically. 
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For L min = 18 we find (£/£)* = 0.56405(5) and j3 c - j3 s = -0.00000067(38). Adding a further correction, we get 
acceptable values of x 2 /DOF already for L min — 7. But also here % 2 /DOF still further decreases and (^nd/L)* and 
p c — (3 S move monotonically with increasing L m i n . For e = 1.6, we obtain the results fand/L)* — 0.56386(10) and 
Pc- P s = -0.00000119(48) for L min = 12. For e = 2 and L min = 12, we get the results (&nd/L)* = 0.56391(8) 
and p c — p s = —0.0000011(46). For the Binder cumulant similar results can be found. We arrive at the final results 
P C (X = 4.5) = 0.6862368(10) and 

(£/£)* = 0.5639(2), Ul = 1.1394(3). (B5) 

The error-bars are chosen such that the results of the different fits are covered. 

A similar analysis for the 0(4) symmetric </> 4 model at A = 12.5 leads to estimates U% = 1.0942(3), £/L = 0.5471(3), 
and P c = 0.909517(2). 



2. Determination of A* 



Next we determine the value of A* where leading corrections to scaling vanish. To this end we study 

Ui(L) = U 4 (L,p f ) (B6) 



= 0.5644 (B7) 



where pf is determined by the equation 



where 0.5644 is the result for (£/£)* of Ref. [6j. In order to compute U 4 we use the first order Taylor expansion (|B4|) 
of £/L and U4 around the simulation point p s . For L = 12 , A = 4.5 we simulate at a number of different P s , to check 
whether this approximation is sufficient for our purpose. In particular we find that for A = 4.5 the difference between 
Ps = 0.686238 and Pf is sufficiently small that contributions oc (P — p s ) 2 can be ignored. Due to scaling, we expect 
that this also holds for all of the lattice sizes that we have simulated. 

First we fit our data obtained at A = 4.5 with a number of different ansaetze 

-0.79 



U A = Ul + clL' ™ , (B8) 



U 4 = Ul + aL- 79 + bL- Cl , (B9) 

and 

Ui = Ul + aL- 79 + bL- ei + cL-' 2 . (B10) 

Also here we fix cj = 0.79; the final results change only little when we replace it with uj — 0.77. In the case of the 
ansatz (|B9[) we set t\ — 1.6 or 2. Finally in ansatz (|B10|) we add two terms with subleading corrections. We have 
fitted using various choices for t\ and 62- 

In our fits we take into account all lattices sizes L > L m i n . In the case of the ansatz (|B8|) we get an acceptable 
X 2 /DOF starting from L m i„ = 22. From this fit we get a = 0.00254(31). Further increasing L m i n , a is monotonically 
increasing; for L min = 30 we obtain a = 0.0037(6). 

Fitting with the ansatz (|B9[) and t\ = 1.6 we obtain an acceptable x 2 /DOF already starting from L m i n = 6. We 
get a = 0.01038(27) for the correction amplitude. Increasing L m i n the correction amplitude remains stable. Using 
instead ei = 2 we get an acceptable x 2 /DOF starting from L min — 7. The corresponding result for the correction 
amplitude is a = 0.00586(20). Increasing L min , the value of a increases up to a = 0.00676(33) for L min = 10. For 
L min = 11 and 12 we get a very similar result. For L min = 12, x 2 /DOF = 14.50/21 and 15.78/21 for e = 2 and 1.6, 
respectively. 

Finally we fit with the ansatz (|B9[) using (ei, 62) = (1.6, 2), (1.6, 1.96) or (1.8, 2). The results of such fits are all in 
the interval 0.005 < a < 0.011. We conclude a — 0.007(4), where the central value and the error-bar are chosen such 
that the results of the different fits are covered. Next we convert this estimate of the correction amplitude at A = 4.5 
into a new estimate of A*. In order to compute the derivative of a with respect to A, we study the differences 



AU 4 (L) = U 4 (L, A = 5) - U 4 (L, A = 4) 
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TABLE V: Fits with the ansatz (fB12| . 0(3) universality class 

L mln c uj x 2 /DOF 

6 -0.0109(2) 0.775(9) 6.44/9 

7 -0.0109(3) 0.777(12) 6.38/8 

8 -0.0111(4) 0.784(16) 5.86/7 



TABLE VI: Fits of AU 4 (L,22, 12.5) with the ansatz (|BT2|l . 0(4) universality class. 





c 


UJ 


x 2 /dof 


6 


-0.00776(14) 


0.777(8) 


6.77/11 


7 


-0.00764(18) 


0.771(10) 


5.81/10 


8 


-0.00753(22) 


0.765(11) 


5.15/9 


9 


-0.00741(28) 


0.759(15) 


4.64/8 



In this difference U% exactly cancels. Furthermore subleading corrections should cancel to a large extend. Therefore 
we fit our data with the ansatz 

AU 4 (L) = cL~ u . (B12) 

Results of such fits with c and uj as free parameters are given in table [Vj Already starting from L m j n = 6 we get 
an acceptable x 2 /DOF. Furthermore, the value obtained for uj is fully consistent with the field theoretic estimates 
uj = 0.782(13) and ui — 0.794(18) obtained by the perturbative expansion in three dimensions fixed and the e- 
expansion, respectively The facts that x 2 /DOF is small and the result for uj is consistent with the field-theoretical 
ones confirm our assumption that already for the lattice sizes that we consider, AU^L) is dominated by the leading 
correction. 

Fitting with u> = 0.79 fixed, to be consistent with the analysis of U4 at A = 4.5 above, we find c = —0.01126(4) and 
X 2 /DOF= 6.0/8 for L m i n — 8. The result for c changes little, when L m i n is varied. In order to check how well the 
derivative of a with respect to A is approximated by the finite difference, we also have fitted Ui(L, A = 5) — U^L, A = 
4.5). Here we find c = —0.00506(4) and x 2 /DOF= 5.4/8 for L min — 8. Also here, the result for c changes little, when 



da \ 1 



Using these results we arrive at 

A* « 4.5 - a(A = 4.5) I ^ j = 4.5 - 0.007(4)/(-2 x 0.00506(4)) w 5.2(4) . (B13) 
We perform a similar analysis in the case of the 0(4) universality class. Here /3/ is given by 

COM/) 



L 



0.547 (B14) 



where 0.547 is the result for (£/£)* of ref. j8|. First we have analyzed the data for U4 at A = 12.5. The analysis 
is done in much the same way as discussed above in detail for the 0(3) universality class. Fixing uj = 0.79 we find 
a = 0.007(5) as amplitude of the leading correction. 
Next we study the difference 

Af7 4 (i> ) A 1 ,A 2 )=C/4(L,A 1 )-C74(i,A 2 ) . (B15) 

We perform fits for Ai = 22, A2 = 12.5 and Ai — 00, A2 = 12.5 using the ansatz (|B12[) with c and cj as free parameters. 
The results for Ai = 22 and Ai = 00 are given in tables IVT1 and IVTll respectively. 

These results can be compared with ui — 0.774(20) and uj = 0.795(30) from the perturbative expansion at three 
dimensions fixed and the e-expansion, respectively [4| . 

Fixing uj = 0.79 we obtain c = -0.00800(2) (with x 2 /DOF=9. 77/12) as amplitude for the differences Ai = 22 and 
A 2 = 12.5 with L mm = 6 Taking data only for L = 6 we get c(Ai = 14, 12.5) = -0.00193(5) c(Ai = 20, 12.5) = 
-0.00696(5), c(Ai = 30,12.5) = -0.01084(5) c(Ai = 32,12.5) = -0.01132(5). It is quite clear from these numbers 
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TABLE VII: Fits of AU 4 (L, oo, 12.5) with the ansatz (|B12[) , 0(4) universality class. 



Lmiu c lu X 2 /DOF 

6 -0.01870(16) 0.787(4) 11.64/7 

7 -0.01849(21) 0.783(5) 9.39/6 

8 -0.01841(26) 0.781(6) 9.05/5 

9 -0.01844(32) 0.782(7) 9.02/4 

10 -0.01846(38) 0.782(8) 9.00/3 

11 -0.01777(45) 0.769(10) 2.03/2 



that a linearization of the correction amplitude as a function of A is not sufficient to compute the estimate of A* . For 
the same reason, we give an asymmetric estimate of the error: 

A* = 20±* 5 ( B16 ) 

This value is larger than A* = 12.5(4.0) that we quote in ref. 8]. However we are quite confident that indeed a A* 
exists for the 0(4) case. Note that in the limit N — > oo for the simple cubic lattice and the given lattice action, no 
A* exists and that leading corrections are minimal in the limit A — > oo [l|. 



3. The magnetic susceptibility and the exponent r\ 

In order to obtain the critical exponent 77, we analyze the behavior of 

X = X(M (B17) 

where, in the 0(3) case (3f is defined by ^(j3f)/L — 0.5644. In the first step of the analysis we eliminate leading 
corrections to scaling. To this end we analyze the ratios 

X( A = 5) _ _„ , „ r - .79 



X(A - 4) 



a(l + cL- u - 79 ) . (B18) 



We obtain a good fit starting from L min = 11. For L min — 11 we obtain a = 0.99172(8), c = —0.0046(6) and 
X 2 /DOF= 3.11/4. Therefore in order to eliminate corrections at A = 4.5 we follow the strategy discussed in section 
IIIIA1I Using U 4 = Ul + 0.007(4)£-°- 79 + ... and U A {\ = 5) - U 4 (X = 4) = -0.01126(4)£-°- 79 ...) Eq. ([29]) reads 



We fit x with the ansaetze 



^^I'-^S - 007 ' 4 '^ (B19) 

X = aL 2 - r i , (B20) 
X = aL 2 - r ' + c , (B21) 

X = aL 2 - ri (l + bL- e ) + c (B22) 

with e = 1.6 or e = 1.8. In the case of the ansatz (|B20[) we obtain very large % 2 /DOF up to L m i n = 32. For 
L mm = 48 we get r\ = 0.0375(1) and x 2 /DOF= 0.46/1. Using the ansatz (|B%|) we get x 2 /DOF« 1 alrea dy for 
L mm = 16; for example, for L mm = 18 we obtain 77 = 0.03767(4) and x 2 /DOF= 10.11/15. Using the ansatz (|B22j) 
with e = 1.6 we get for L min = 10 the results 77 = 0.03791(7) and x 2 /DOF= 14.55/22. and for e = 1.8 and L min = 8 
we get i] = 0.03780(3) and x 2 /DOF= 18.74/24. We redo these fits for \ without correction to check the effect of 
the uncertainty of A*. We find that the estimates of rj change by about 0.0001. Taking into account only fits with 
ansaetze that include the analytic background, we arrive at 

77 = 0.0378(3) . (B23) 

In the case of the 0(4) universality class, performing a similar analysis we obtain 

77 = 0.0360(3) . (B24) 
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4. The exponent v 



We estimate the exponent v from the behavior of the slope of C/4 and £/L at f3 c : 




aL^^ + cL'" + ...)• 



(B25) 



Since we did not plan to compute the exponent v from the beginning, we did not compute the second derivatives of 
C/4 and £/L with respect to f3. Hence we can not compute the slope at fixed values of C/4 or £,/L. At A = 4.5 we 
performed MC simulation very close to our final value of (3 C . Therefore it is sufficient to have a rather rough estimate 
of the second derivatives of C/4 and £/L in order to compute the first derivatives of C/4 and £/L at (3 C starting from 
the first derivatives of C/4 and £/L at j3 s that we have computed in our simulations. To this end, we simulated for 
L = 12 at a number of different (3 values. Using these data we compute the second derivatives of C/4 and £/L with 
respect to /3 by finite differences. The second derivatives are then estimated by R"(L) = R" '(12)(L /12) 2 / v . Notice 
that our estimate of j3 c = 0.6862368(10) is very close to the simulation point j3 s = 0.686238. We analyze the resulting 
data by fitting with various ansaetze that are derived from Eq. (|B25|) . We arrive at v = 0.7118(7) from the analysis 
of the slope of of £/L and v = 0.7114(11) from that of C/4. The error bars take also into account the uncertainty of 
A*. As our final estimate we quote 
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